

/*
plot qe simulation output

output:
- f5c_*: model experiment figures

*/


* plotting options
set scheme lean_uncluttered


* command line option
* set to 1 to plot all robustness experiments
if "`1'"=="" {
	local plot_model_exp 0
}
else {
	local plot_model_exp `1'
}



********************************************************************************
********************************************************************************
* simple program to plot model experiments
cap program drop plot_qe_experiments 
program define plot_qe_experiments

* parameter arguments
syntax, baseline_val(real) param_name(varname) ///
	xaxis_title(string) ///
	[plot_sd(real 0) experiment_name(string) ///
	legend_off(real 0) legend_pos(real 2)]

if "`experiment_name'" == "" {
	local experiment_name "qe_experiment"
}
if `legend_off' == 0 {
	local legend_opts order(1 "Baseline" 2 "Risky-only" 3 "Safe-only" ) ///
		ring(0) position(`legend_pos') rows(3) ///
		size(small) rowgap(0) colgap(0)
}
else {
	local legend_opts off
}
if `plot_sd' == 0 {
* inflation (long-run level)
twoway line pi_oo pi_oo_riskyonly pi_oo_safeonly `param_name', ///
	xline(`baseline_val', lcolor(black) lpattern(dot) ) ///
	/// xtitle("`xaxis_title'") ///
	/// ytitle("Inflation Response" "(Long-Run)") ///
	xtitle("") ytitle("") ///
	legend(`legend_opts' ) ///
	name("`experiment_name'_pi", replace) nodraw
* output (long-run level)
twoway line x_oo x_oo_riskyonly x_oo_safeonly `param_name', ///
	xline(`baseline_val', lcolor(black) lpattern(dot) ) ///
	/// xtitle("`xaxis_title'") ///
	/// ytitle("Output Response" "(Long-Run)") ///
	xtitle("") ytitle("") ///
	legend(`legend_opts' ) ///
	name("`experiment_name'_x", replace) nodraw
}
else {
* inflation (long-run sd)
twoway line sd_pi_oo sd_pi_oo_riskyonly sd_pi_oo_safeonly `param_name', ///
	xline(`baseline_val', lcolor(black) lpattern(dot) ) ///
	/// xtitle("`xaxis_title'") ///
	/// ytitle("Inflation Volatility" "(Long-Run)") ///
	xtitle("") ytitle("") ///
	legend(`legend_opts' ) ///
	name("`experiment_name'_pi", replace) nodraw
* output (long-run sd)
twoway line sd_x_oo sd_x_oo_riskyonly sd_x_oo_safeonly `param_name', ///
	xline(`baseline_val', lcolor(black) lpattern(dot) ) ///
	/// xtitle("`xaxis_title'") ///
	/// ytitle("Output Volatility" "(Long-Run)") ///
	xtitle("") ytitle("") ///
	legend(`legend_opts' ) ///
	name("`experiment_name'_x", replace) nodraw
}

end



********************************************************************************
********************************************************************************
* baseline experiments: QE1 comparison
n di "plotting QE1 baseline experiment"

use ../data/fred/_fred_treasury_rates.dta, clear
drop DGS*MO
sort date
gen _T = _n
tsset _T
ds date _T, not
foreach var of varlist `r(varlist)' {
	replace `var' = 100*`var'
	gen D_`var' = F.`var' - L.`var'
}
keep date D_*
* QE1 dates
keep if date==mdy(11,25,2008) | ///
	date==mdy(12,1,2008) | ///
	date==mdy(12,16,2008) | ///
	date==mdy(1,28,2009) | ///
	date==mdy(3,18,2009)
drop date
collapse (sum) *
gen _i = 1
reshape long D_DGS , i(_i) j(maturity)
drop _i

rename maturity tau

* merge model-implied response
merge 1:1 tau using ../data/model/experiments_QE1_yc.dta, nogen
* basis points
replace y_tau = 100*y_tau
sort tau

twoway (scatter D_DGS tau, m(Oh)) (line y_tau tau) , ///
	xtitle("Maturity") ytitle("Response (b.p.)") ///
	legend(region(fcolor(none)) order(1 "Data" 2 "Model") ///
		ring(0) position(2) rows(2) ) ///
	yscale(r(-115 0) titlegap(-5)) ylabel(0[-25]-100) ///
	name("qe1_fit", replace) 
local fname f5c_qe1_qt_fit_irf_p1
graph export ../figs/`fname'.pdf, replace
graph export ../figs/`fname'.eps, replace


* QE1 IRFs
use ../data/model/experiments_QE1_irf.dta, clear

* IRF (output and inflation)
twoway line x pi t if t<=5, ///
	xtitle("Time (Years)") ///
	ytitle("Response (ppt.)") ///
	legend(region(fcolor(none)) order(1 "Output" 2 "Inflation" ) ///
		ring(0) position(2) rows(2) ) ///
	yscale(titlegap(0)) ///
	name("qe1_irf", replace) 
local fname f5c_qe1_qt_fit_irf_p2
graph export ../figs/`fname'.pdf, replace
graph export ../figs/`fname'.eps, replace


* QT IRF
use ../data/model/experiments_QT_passive_med_irf.dta, clear

* IRF (output and inflation)
twoway line x pi t if t<=5, ///
	xtitle("Time (Years)") ///
	ytitle("Response (ppt.)") ///
	ysc(titlegap(20)) ///
	legend(region(fcolor(none)) order(1 "Output" 2 "Inflation" 3 " ") ///
		ring(0) position(5) cols(1) ) ///
	yscale(titlegap(-15)) ///
	name("qt_irf", replace) 
local fname f5c_qe1_qt_fit_irf_p3
graph export ../figs/`fname'.pdf, replace
graph export ../figs/`fname'.eps, replace



* combine with correct size and export
local fname f5c_qe1_qt_fit_irf
graph combine qe1_fit qe1_irf qt_irf, ///
	rows(3) imargin(small) xsize(6) ysize(10) 
*scale(0.6)

* add titles
gr_edit .plotregion1.graph1.title.text.Arrpush "Panel A: Yield Response to QE1"
gr_edit .plotregion1.graph2.title.text.Arrpush "Panel B: Macroeconomic Response to QE1"
gr_edit .plotregion1.graph3.title.text.Arrpush "Panel C: Macroeconomic Response to QT"
graph export ../figs/`fname'.pdf, replace
graph export ../figs/`fname'.eps, replace




********************************************************************************
********************************************************************************
* compare IRF of QE1 and QT shocks
n di "plotting QE1/QT comparison"

use ../data/model/experiments_QE1_irf.dta, clear
keep bQE t
rename bQE bQE1

merge 1:1 t using ../data/model/experiments_QT_passive_med_irf.dta, ///
	nogen keepusing(bQE bQE_p)
rename bQE bQT_a
rename bQE_p bQT_p

gen bQT = 0.25*bQT_a + 0.75*bQT_p


* normalize
gen bQE_norm = bQE1 / bQE1[1]
gen bQT_norm = bQT / (-bQE1[1])
gen bQT_a_norm = 0.25*bQT_a / (-bQE1[1])
gen bQT_p_norm = 0.75*bQT_p / (-bQE1[1])

twoway line bQE_norm bQT_norm bQT_a_norm bQT_p_norm t, ///
	legend( ///
		order(1 "QE1" 2 "QT" 3 "QT: active" 4 "QT: passive" ) ///
		ring(0) position(2) rows(4) ///
	) xtitle("Time (Years)") ytitle("Magnitude (Relative to QE1)") ///
	xsize(8) ysize(5) name("qt_qe1_comparison", replace) 

* save
local fname f5c_qt_qe1_comparison
graph export ../figs/`fname'.pdf, replace
graph export ../figs/`fname'.eps, replace




************************************************************
************************************************************
if `plot_model_exp'==1 {
n di "plotting QE robustness"
* plot model experiments

**** main experiments ****
* risk aversion
local experiment_name qe_exp_risk_aversion
use ../data/model/`experiment_name'.dta, clear
* normalize
replace a = a / 0.1
plot_qe_experiments , baseline_val(1.0) plot_sd(0) param_name(a) ///
	xaxis_title("Risk Aversion") legend_pos(11) ///
	experiment_name("`experiment_name'")

* mean reversion
local experiment_name qe_exp_mean_reversion
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(0.2) plot_sd(0) param_name(kappa_QE) ///
	xaxis_title("Mean Reversion") legend_off(1) ///
	experiment_name("`experiment_name'")

* payoff risk
local experiment_name qe_exp_payoff_risk
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(1.115727778) plot_sd(0) param_name(s_d) ///
	xaxis_title("Risky Payoff Volatility") legend_pos(11) legend_off(1) ///
	experiment_name("`experiment_name'")

* QE volatility
local experiment_name qe_exp_policy_volatility
use ../data/model/`experiment_name'.dta, clear
* normalize
replace sigma_QE = sigma_QE / 5.122532344
plot_qe_experiments , baseline_val(0.0) plot_sd(1) param_name(sigma_QE) ///
	xaxis_title("QE Volatility (Relative to Demand Shocks)") legend_pos(11) legend_off(1) ///
	experiment_name("`experiment_name'")

* QE maturity duration
local experiment_name qe_exp_maturity_dur
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(0.35) plot_sd(0) param_name(theta1_QE) ///
	xaxis_title("QE Maturity Duration") legend_pos(7) legend_off(1) ///
	experiment_name("`experiment_name'")

* effective borrowing shape
local experiment_name qe_exp_effrate_shape
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(2.0) plot_sd(0) param_name(eta1_til) ///
	xaxis_title("Eff. Borrowing Maturity Weights") legend_pos(2) legend_off(1) ///
	experiment_name("`experiment_name'")


* taylor rule: output
local experiment_name qe_exp_mp_x
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(0.0) plot_sd(0) param_name(p_ix) ///
	xaxis_title("Taylor Rule: Output") legend_pos(2) legend_off(1) ///
	experiment_name("`experiment_name'")

* taylor rule: inflation
local experiment_name qe_exp_mp_pi
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(3.096074726) plot_sd(0) param_name(p_ipi) ///
	xaxis_title("Taylor Rule: Inflation") legend_pos(2) legend_off(1) ///
	experiment_name("`experiment_name'")

* taylor rule: inertia
local experiment_name qe_exp_mp_reversion
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(1.081595053) plot_sd(0) param_name(k_i) ///
	xaxis_title("Taylor Rule: Inertia") legend_pos(2) legend_off(1) ///
	experiment_name("`experiment_name'")


**** additional experiments ****
* price rigidity
local experiment_name qe_exp_pc_rigidity
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(0.705042668) plot_sd(0) param_name(delta_x) ///
	xaxis_title("Price Rigidity") legend_pos(2) legend_off(1) ///
	experiment_name("`experiment_name'")

* dividend process: output
local experiment_name qe_exp_dividend_x
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(0.39270544) plot_sd(0) param_name(p_dx) ///
	xaxis_title("Dividend Process: Output") legend_pos(7) /// legend_off(1) ///
	experiment_name("`experiment_name'")

* dividend process: inflation
local experiment_name qe_exp_dividend_pi
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(0.0) plot_sd(0) param_name(p_dpi) ///
	xaxis_title("Dividend Process: Inflation") legend_pos(2) legend_off(1) ///
	experiment_name("`experiment_name'")

* dividend process: inertia
local experiment_name qe_exp_dividend_reversion
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(0.879102025) plot_sd(0) param_name(k_d) ///
	xaxis_title("Dividend Process: Inertia") legend_pos(2) legend_off(1) ///
	experiment_name("`experiment_name'")

* effective borrowing weight
local experiment_name qe_exp_effrate_weight
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(0.0) plot_sd(0) param_name(n0) ///
	xaxis_title("Eff. Borrowing Risky/Safe Weights") legend_pos(2) legend_off(1) ///
	experiment_name("`experiment_name'")

* QE active passive
local experiment_name qe_exp_active_passive
use ../data/model/`experiment_name'.dta, clear
plot_qe_experiments , baseline_val(0.0) plot_sd(0) param_name(active_passive_frac) ///
	xaxis_title("QE Active/Passive Fraction") legend_pos(3) legend_off(1) ///
	experiment_name("`experiment_name'")


******************************
* combined figs

* main policy experiments
* inflation
local fname f5c_qe_model_exercises_infl
graph combine ///
	qe_exp_risk_aversion_pi qe_exp_payoff_risk_pi ///
	qe_exp_mean_reversion_pi qe_exp_maturity_dur_pi ///
	qe_exp_effrate_shape_pi qe_exp_effrate_weight_pi ///
	qe_exp_mp_pi_pi qe_exp_mp_x_pi ///
	qe_exp_mp_reversion_pi qe_exp_policy_volatility_pi, ///
	rows(5) imargin(tiny) xsize(12) ysize(16) name(`fname', replace)
* add titles
gr_edit .plotregion1.graph1.title.text.Arrpush "Panel A: Risk Appetite"
gr_edit .plotregion1.graph2.title.text.Arrpush "Panel B: Payoff Risk"
gr_edit .plotregion1.graph3.title.text.Arrpush "Panel C: QE Unwinding"
gr_edit .plotregion1.graph4.title.text.Arrpush "Panel D: QE Asset Duration"
gr_edit .plotregion1.graph5.title.text.Arrpush "Panel E: Eff. Wgt. Maturity"
gr_edit .plotregion1.graph6.title.text.Arrpush "Panel F: Eff. Wgt. Riskiness"
gr_edit .plotregion1.graph7.title.text.Arrpush "Panel G: Mon. Pol. Inflation"
gr_edit .plotregion1.graph8.title.text.Arrpush "Panel H: Mon. Pol. Output"
gr_edit .plotregion1.graph9.title.text.Arrpush "Panel I: Mon. Pol. Inertia"
gr_edit .plotregion1.graph10.title.text.Arrpush "Panel J: QE Volatility"
graph export ../figs/`fname'.pdf, replace
graph export ../figs/`fname'.eps, replace

* output
local fname f5c_qe_model_exercises_output
graph combine ///
	qe_exp_risk_aversion_x qe_exp_payoff_risk_x ///
	qe_exp_mean_reversion_x qe_exp_maturity_dur_x ///
	qe_exp_effrate_shape_x qe_exp_effrate_weight_x ///
	qe_exp_mp_pi_x qe_exp_mp_x_x ///
	qe_exp_mp_reversion_x qe_exp_policy_volatility_x, ///
	rows(5) imargin(tiny) xsize(12) ysize(16) name(`fname', replace)
* add titles
gr_edit .plotregion1.graph1.title.text.Arrpush "Panel A: Risk Appetite"
gr_edit .plotregion1.graph2.title.text.Arrpush "Panel B: Payoff Risk"
gr_edit .plotregion1.graph3.title.text.Arrpush "Panel C: QE Unwinding"
gr_edit .plotregion1.graph4.title.text.Arrpush "Panel D: QE Asset Duration"
gr_edit .plotregion1.graph5.title.text.Arrpush "Panel E: Eff. Wgt. Maturity"
gr_edit .plotregion1.graph6.title.text.Arrpush "Panel F: Eff. Wgt. Riskiness"
gr_edit .plotregion1.graph7.title.text.Arrpush "Panel G: Mon. Pol. Inflation"
gr_edit .plotregion1.graph8.title.text.Arrpush "Panel H: Mon. Pol. Output"
gr_edit .plotregion1.graph9.title.text.Arrpush "Panel I: Mon. Pol. Inertia"
gr_edit .plotregion1.graph10.title.text.Arrpush "Panel J: QE Volatility"
graph export ../figs/`fname'.pdf, replace
graph export ../figs/`fname'.eps, replace



* additional robustness
* inflation
local fname f5c_qe_model_robustness_infl
graph combine ///
	qe_exp_dividend_x_pi qe_exp_dividend_pi_pi ///
	qe_exp_dividend_reversion_pi qe_exp_active_passive_pi, ///
	rows(2) imargin(tiny) xsize(12) ysize(8) name(`fname', replace)
* add titles
gr_edit .plotregion1.graph1.title.text.Arrpush "Panel A: Payoff Risk: Output"
gr_edit .plotregion1.graph2.title.text.Arrpush "Panel B: Payoff Risk: Inflation"
gr_edit .plotregion1.graph3.title.text.Arrpush "Panel C: Payoff Risk: Inertia"
gr_edit .plotregion1.graph4.title.text.Arrpush "Panel D: QE Active/Passive"
graph export ../figs/`fname'.pdf, replace
graph export ../figs/`fname'.eps, replace

* output
local fname f5c_qe_model_robustness_output
graph combine ///
	qe_exp_dividend_x_x qe_exp_dividend_pi_x ///
	qe_exp_dividend_reversion_x qe_exp_active_passive_x, ///
	rows(2) imargin(tiny) xsize(12) ysize(8) name(`fname', replace)
* add titles
gr_edit .plotregion1.graph1.title.text.Arrpush "Panel A: Payoff Risk: Output"
gr_edit .plotregion1.graph2.title.text.Arrpush "Panel B: Payoff Risk: Inflation"
gr_edit .plotregion1.graph3.title.text.Arrpush "Panel C: Payoff Risk: Inertia"
gr_edit .plotregion1.graph4.title.text.Arrpush "Panel D: QE Active/Passive"
graph export ../figs/`fname'.pdf, replace
graph export ../figs/`fname'.eps, replace



}







